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Abstract 

We study cosmic-ray acceleration in young Type la Supernova Remnants (SNRs) by means 
of test-particle diffusive shock acceleration theory and 1-D hydrodynamical simulations of 
their evolution. In addition to acceleration at the forward shock, we explore the particle 
acceleration at the reverse shock in the presence of a possible substantial magnetic field, 
and consequently the impact of this acceleration on the particle spectra in the remnant. 
We investigate the time evolution of the spectra for various time-dependent profiles of the 
magnetic field in the shocked region of the remnant. We test a possible influence on particle 
spectra of the Alfvenic drift of scattering centers in the precursor regions of the shocks. 
In addition, we study the radiation spectra and morphology in a broad band from radio 
to gamma-rays. It is demonstrated that the reverse shock contribution to the cosmic-ray 
particle population of young Type la SNRs may be significant, modifying the spatial distri- 
bution of particles and noticeably affecting the volume-integrated particle spectra in young 
SNRs. In particular spectral structures may arise in test-particle calculations that are often 
discussed as signatures of non-linear cosmic-ray modification of shocks. Therefore, the spec- 
trum and morphology of emission, and their time evolution, differ from pure forward-shock 
solutions. 

Keywords: Supernova Remnants; cosmic rays; cosmic-ray acceleration; hydrodynamics; 
forward and reverse shocks 



1. Introduction 

The energy density in local cosmic rays 
(CRs), when extrapolated to the whole 
Galaxy, implies the existence of powerful 
accelerators inside the Galaxy. Supernova 
Remnants (SNRs) have long been thought 
to be the main candidates on account of 
energetics Later, diffusive shock ac- 
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celeration (DSA) of charged particles at 

the SNR shocks was considered i a a 

and is now thought to be the most likely 
mechanism that can accelerate CRs up to 
10^^ eV, where the cosmic-ray spectrum 
shows a break known as the 'knee' If 
they are very efficiently accelerated, CRs 
can produce a feedback on the plasma flow 
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in SNRs through their pressure, and can 
modify the shock structure by decelerating 
the incoming plasma flow when streaming 
through it in the precursor region (see [3] 
for a review and references therein). Ob- 
servationally, it is unclear whether or not 
the CR-acceleration efficiency is sufficiently 
high for a significant feedback jl, ^ , but if 
so, SNR shocks may become modified, and 
the particle spectra may deviate from the 
classical power-law slope of s = 2 found 
in test-particle calculations. The spectra 
produced by modified shocks show concave 
structure with softer (s > 2) low-energy and 
harder (s < 2) high-energy parts. 

The CR spectrum observed at Earth sug- 
gests that sources produce softer spectra 
than predicted by non-linear DSA (NDSA) 
theory [lO| . If the CR sources produce very 
hard spectra, the energy dependence of the 
diffusion coefficient in the Galaxy needs to 
be much stronger than is allowed by the ob- 



which increases MHD turbulence 14, 25, 26 



served CR anisotropy [Ul]. It should be 



noted, though, that the CR source spec- 
trum is that of the CR escaping from their 
sources, the time integral of which could 
be considerably softer than the spectrum 
of particles contained in the sources at a 
given time 
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On the other hand, re- 
cent 7-ray observations of some shell-type 
SNRs tend to show soft spectra (F > 2): 
J1713.7-3946jl3,l3, RX J0852.0-4622 
li, RCW 86 [37 SN 1006 [13, Cas A 



M 119|, IC 443 H, Mi and other (22|, |23|, 
12*^ . This means that the particle spectra, 
assuming the radiation was produced by 
hadronic scenario, is also softer than NDSA 
predictions. 

It is known that CR particles are scat- 
tered by magnetohydrodynamic (MHD) 
waves. These waves, alongside amplified 
magnetic fields (MFs), may be generated 
by instabilities invoked by CR streaming. 



27|, or directly by pure MHD turbulence in- 



duced by propagation of coUisionless shock 



in media 28 



29, 30 



Gyro-resonant instabilities favor the gen- 
eration of MHD waves that propagate 



with Alfven velocity [31|, whereas non- 
resonant instabilities can initially favor non- 
oscillating and non-propagating modes. If 
the former dominate, the scattering centers 
move relative to the plasma, and the in- 
clusion of this so-called Alfvenic drift may 
soften the spectra so that they could match 



experimental data [32|, |33|, l3^, |35|, l36(| . Alter- 



natively, the bulk of Galactic CRs may be 
produced in a low-efficiency or test-particle 
mode transferring only < 10% of the SNR 



energy to CRs [3J|. This number is roughly 
that needed to explain the origin of CRs by 
acceleration in Supernova Remnants (37| - 

Given the complexities faced by NDSA 
theory, in this paper we propose to inves- 
tigate the diversity of particle spectra that 
can be produced by young SNR shocks in 
a test-particle mode, concentrating for sim- 
plicity on Type la SNRs in this paper. 
We consider acceleration by both forward 
and reverse shocks in evolving spherically- 
symmetric systems. Many pieces of obser- 
vational evidence suggest that the reverse 
shock (RS) may also accelerate particles to 
high energies. This includes the detection of 
non-thermal X-rays from the RS of Cas A 



38 
41 



39 



40 , and RCW 86 



IE 0102.2-7219 
as well as radio emission from the RS 
CR acceleration at 



of Kepler's SNR [42 



the RS is also supported by theory |43|, |44 



In any case, high-energy particles produced 
by the^FS may be re-accelerated at the 
RS. 
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The preponderance of observa- 
tional and experimental evidence, combined 
with earlier theoretical implications, sug- 
gests that it is essential to at least consider 



the effects of botli sliocks in accelerating 
particles. The total volume-integrated spec- 
tra in such will deviate from expec- 
tations based on plane-parallel shocks and 
steady-state systems employing the FS only. 

The MF in the shocked interaction re- 
gion is essential to determining the parti- 
cle spectra, but is extremely hard to deduce 
observationally. This has led to the param- 
eterization of the field in different models, 
which generally involves tying the field to 
dynamically or kinematically determinable 
parameters. Our intent is not to study the 
amplification of the field itself, which is be- 
yond the scope of this paper, but to investi- 
gate the particle spectra that are produced 
under widely used assumptions of the field 
in the literature. We also explore how the 
inclusion of Alfvenic drift affects the total 
particle spectra in SNR. Given the obser- 
vational data on MF amplification 46l . |47 
as well as theoretical M, UM and numerical 



48L \2& results, we consider two cases: 



moderate (75/iG) and a strongly amplified 
(300/iG) MF. Besides calculating the to- 
tal particle spectra, we produce wide-band 
spectral-energy distributions (SEDs), and 
provide surface brightness profiles in the ra- 
dio. X-ray and 7-ray bands. 

This work intends to improve upon the 
present state of the field in many ways 

(1) We take into account the free expan- 
sion phase of the expanding SNR, when the 
shock velocities are the highest, whereas 
most (though certainly not all) analyses 
have totally neglected the initial phase, 
and started directly with the Sedov phase. 

(2) Our technique to solve the transport 
equation, and the corresponding high res- 
olution spherically-symmetric simulations 
used, guarantee a very high level of ac- 
curacy and therefore an accurate render- 
ing of the particle spectra for the chosen 



MF and diffusion coefficient. (3) We have 
taken into account the acceleration of CRs 
at the reverse shock, as suggested by re- 
cent observations. We have shown what ob- 
servable features in SNRs we may expect 
if a non-negligible MF permits CR accel- 
eration at the RS to high energies. This 
has not been explored in such detail be- 
fore, including spherical symmetry and fully 
time-dependent cosmic-ray transport. Ear- 
lier studies |43|] investigated CR acceleration 
at the RS using an approximate, analyti- 
cal model [iij and assumed there is no MF 
amplification at the RS. In particular, they 
focused on non-linear effects arising from a 
high injection efficiency in low MF. Other 
authors [4J] concentrated on the study of 
a specific object, which most probably is 
a core-collapse SNR, and have not made 
any general conclusions regarding type-la 
SNRs. (4) Our results therefore present a 
more complete study of the particle spec- 
tra and high-energy emission that is to be 
expected from type-la SNe with moderate 
acceleration efficiency. 



2. Cosmic-Ray Acceleration 



We consider a test-particle approach to 
diffusive shock acceleration of CRs in SNRs 
0, i, B, H- It has been recently jij shown 
that a test-particle description is applica- 
ble if the CR pressure at the shock is less 
than 10% of the shock ram pressure. There- 
fore, if one limits the amount of energy con- 
tained in CRs and the CR pressure at the 
shock, the acceleration can be treated by 
independently solving the cosmic-ray trans- 
port equation and the hydrodynamic equa- 
tions of SNR evolution. The cosmic-ray 
transport equation is a diffusion-advection 



equation in both space and momentum 3l| : 
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where is the differential number density 
of cosmic rays, is the spatial diffusion co- 
efficient, V is the advective velocity given by 
a 1-D hydrodynamical simulation, p are the 
energy losses, and Q is the source term rep- 
resenting the injection of the thermal parti- 
cles into the acceleration process given as 

Q = Vinu\Vsh-~Vu\6{r-Rsh)S{p-pinj), (2) 

where rji is the injection efficiency parame- 
ter, Uu is the number density of plasma in 
the shock upstream region, Vsh is the shock 
speed, Vu is the plasma velocity in the shock 
upstream region, r is the distance from the 
SNR center, Rgh is the radius of the shock, 
p is the particle momentum, and pmj is the 
momentum of the injected particles. 

We assume the thermal leakage injection 
model 



50 



where only the particles with 
momentum pinj > ippth can be accelerated, 
with ip being a multiple of the particle ther- 
mal momentum, pth- The efficiency of in- 
jection is determined as 



30F 



(3) 



where a is the shock compression ratio. Al- 
though a and pth are obtained from hydro 
simulations, ip remains a free parameter in 
the model. In our calculations we assume 
ip ~ 4.45, keeping rji c:^ 5 x 10~^ sufficiently 
low to stay within the framework of the test- 
particle approximation. We note that since 
we are not interested here in the absolute ra- 
tio of electron to proton spectra, and lacking 
secure knowledge of the details of electron 



injection, as well as the time-scale of ther- 
malization between electrons and protons, 
we assume that electrons and protons are 
injected equally and are at the same tem- 
perature, i.e. the injected electron to pro- 



ton ratio, K, 



e/p,i 



1. Thus, Ke/p obtained 



in our simulations is purely a result of the 
electron-to-proton mass ratio and corrected 
for standard ISM abundances. 

We impose a free-escape spatial bound- 
ary upstream of the FS that should account 
for escape of the highest-energy particles 
from the system. We assume that all par- 
ticles crossing the boundary leave the sys- 
tem, therefore we set the CR number den- 
sity, N, to zero at the free-escape boundary, 
at i?esc = '2'Rfs, where Rps is the FS radius. 

We also assume that the highest en- 
ergy particles generate MHD waves in the 
shock upstream regioii, and thus may am- 
plify the MF H 0, 0, [ij- Being scat- 
tered by MHD waves, the CR particles un- 
dergo Bohm diffusion everywhere inside the 
escape boundary [5l|. Also, the advec- 
tive velocity in Eq. [T], given the presence 
of amplified MF, may include the Alfvenic 
drift of scattering centers and thus may 
not coincide with the plasma flow velocity 



32, 35, 34, 33 



In the current paper we assume spher- 
ical symmetry and do not consider the 
momentum-diffusion term representing 
stochastic (second order Fermi) accelera- 
tion. 

The main difficulty for the numerical so- 
lution of Eq. [T] is that the diffusion coeffi- 
cient Dr is strongly dependent on the parti- 
cle momentum {Dj. = Dr{p) ~ p, in case 
of Bohm diffusion) and therefore spans a 
wide range of magnitudes, ln(pjnax/Pmj) ~ 
25. Formally, the spatial grid must re- 
solve the smallest diffusion length defined 
as Dr{pinj)/Vsh, otherwise the modeled ac- 



celeration is artificially slow. The lowest- 
energy particles have a very small diffusion 
length, and thus a uniform grid must con- 
tain millions (or even billions) of cells which 
is computationally impossible. 

Several approaches are known to date to 
overcome this numerical difficulty. The so- 
called normalized grid is used in 52| , where 
a substitution of variables is performed and 
the spatial coordinates are normalized to 
the diffusion scale. Another approach pre- 
sented by [53| consists of introducing a co- 
moving frame, in which the shock is station- 
ary, coupled with a mesh refinement - a fine 
discretization is used only near the shock 
region, where the low-energy particles are 
injected. 

Here we use the second approach. We 
normalize the spatial coordinate in Eq. [T] 
to the shock radius, introducing coordinate 
X = r/Rsh, and then substitute the coordi- 
nate X with a new coordinate, x*, for which 
we use a uniform grid when solving Eq. [T] 



{x-1) 



with Jacobian 



dx 

doc ^ 



3 (x^ 



(4) 



(5) 



These simple transformations allow us i) 
to place the shock at the center of a known 
cell at any given time and ii) to resolve the 
shock vicinity for low-energy particles with 
only a few hundred bins in the spatial coor- 
dinate, x^. We use parametrizations for the 
radial dependence of the MP as described 
in Section HI and we obtain the hydrody- 
namical parameters (the plasma flow veloc- 
ity, density and the shock speed) needed to 
solve Eq. [T] from 1-D hydrodynamic simu- 
lations of the SNR evolution described in 
Section [31 



Then, Eq. [T] is solved in spherically- 
symmetric geometry using implicit finite- 
difference methods implemented in the FiPy 
54( 1 library modules. We make separate 



runs for the FS and the RS, using the ap- 
propriate coordinate transformations to ob- 
tain the highest resolution where injection 
occurs. Then we sum the results because 
they are linear and independent. Some 
of the highest-energy particles may reach 
the other shock and be re-accelerated there. 
However, the energy of these particles is al- 
ready so high, and their number so low, that 
the poor resolution at the other shock does 
not visibly affect the results. Finally, we 
check the CR-to-ram pressure ratio as well 
as energy containment in CRs a posteriori. 
In our calculations these values are never 
allowed to exceed 10%. 

3. Hydrodynamics 

To study the evolution of the SNR 
shock waves over time, we have performed 
spherically-symmetric hydrodynamic sim- 
ulations using the VH-1 code 55[, a 3- 
dimensional hydrodynamic code based on 
the Piecewise Parabolic Method The 
simulations are based on an expanding grid 
57, 58 , which tracks the motion of the 



outer shock and expands along with it. 
Thus, a high resolution is maintained right 
from the start of the simulations, which is 
essential considering that the remnant size 
increases by about 6 orders of magnitude 
during the simulation. The expanding grid 
also has another advantage - the position of 
the forward shock is almost stationary on 
the grid over most of the evolution, which 
is useful in computing the particle spec- 
trum. Furthermore, the grid expansion is 
adjusted such that, as far as possible, the 
shocked interaction region between the in- 
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ner and outer shocks occupies most of the 
grid. Although coohng of the material is in- 
corporated into the simulations via a cool- 
ing function, the velocities and densities are 
such that cooling is not effective and does 
not play a role. 

In order to study the evolution of young 
SNe in the ambient medium, we need at 
minimum two parameters - a description of 
the density profile of the material ejected 
in the explosion, and that of the medium 
into which the SN is expanding. In this 
first paper, we consider the evolution of a 
Type la SN, which are thought to arise from 
the thermonuclear deflagration of a white 
dwarf. Since the progenitor star is a low- 
mass star and (presumably) does not un- 
dergo significant mass-loss, the medium sur- 
rounding it is assumed to be a constant- 
density interstellar medium. 

The ejecta density structure is more 
uncertain. However, by comparing 
spherically-symmetric models ^_of Type la 
SN explosions, it was found 
ejecta structure of a Type 



59 



that the 
la SN after 
explosion can best be represented by an 
exponential ejecta density profile. This is 
the profile that we have used in our simu- 
lations. Since an exponential introduces an 
additional dimensional parameter to satisfy 
a non-dimensional exponent, the resulting 
SN evolution is no longer self-similar. The 
SN is homologously expanding, therefore 
the ejecta velocity simply increases linearly 
with radius. 

The initial conditions then depend on 
three parameters: (i) The energy in the 
ejected material, which we take as the 
canonical energy of a SN explosion, 10^^ erg, 

(ii) the mass of the ejecta, which we take to 
be the Chandrashekhar mass, 1.4 Mq, and 

(iii) the density of the surrounding medium, 
which we take to be 9.36 xlO^^^ g cm~'^, or 
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Figure 1: Hydrodynamical parameters of the SNR 
at different ages. The x-axis is scaled to the radius 
of the forward shock. 



a number density of about 0.4 for a gas with 
90% H, 10% He, and a trace of other mate- 
rials. 

The supersonic expansion of the SN 
ejecta into the ambient medium gives rise 
to a forward shock expanding into the 
medium, and a reverse shock expanding 
back into the ejecta in a Lagrangian sense. 
The two are separated by a contact discon- 
tinuity, which separates the shocked ejecta 
from the shocked ambient medium. Al- 
though not captured in our spherically sym- 
metric simulations, in mult i- dimensions the 
decelerating contact discontinuity is always 
unstable to Rayleigh- Taylor (R-T) instabil- 
ities, leading to R-T "fingers" of shocked 
ejecta expanding into the shocked ambient 
medium. 

The simulation commences with a grid of 
size 6 xlO~^ pc. By the end of the sim- 
ulation, which is carried on for about 1000 
years, it has grown to more than 10 pc. The 
hydrodynamic parameters at the age of 100, 
400, and 1000 years are shown in Fig[T] The 
outer shock, inner shock and contact discon- 
tinuity are all clearly delineated as sharp 
discontinuities in the density panel (top). 
After a few hundred years the shocked re- 
gion occupies more than 90% of the grid. 
This resolution is essential for computing 
the velocity differential across the shock, 
required in the solution of the cosmic-ray 
transport equation. Note that, while the 
velocity varies smoothly across the contact 
discontinuity (CD), and the pressure shows 
a slight change in slope, there is a dra- 
matic change in density across the CD. The 
shocked ejecta reach a maximum density on 
the inner side of the contact discontinuity, 
while the shocked ambient medium reaches 
a minimum on the outer side. The temper- 
ature therefore reaches a maximum at the 
CD for the shocked ambient medium. This 



distinguishes it from the structure obtained 
for a power-law profile expanding into the 
ambient medium, as would be more appro- 
priate for a Type II or core-collapse SN [59 . 

Finally, the numerical solution of the hy- 
drodynamic equations is mapped onto the 
spatial coordinate x* (cf. EqH]) in which we 
have written the particle transport equation 
(Eq. [1]). The shocks, which are invariably 
smeared out in the diffusive hydrodynamic 
calculation, are re-sharpened by interpola- 
tion toward a point-like jump in density, 
flow velocity, and pressure. This is neces- 
sary to maintain a realistically fast acceler- 
ation at GeV-to-TeV energies. 

4. Magnetic Field and Models Setup 



4.I. Magnetic Field Amplification at SNR 
shocks 

The strength of the magnetic field in the 
SNR is one of the crucial parameters for 
CR acceleration, particularly so in the shock 
regions. Since we assume Bohm diffusion 
here, the MF determines the mean free path 
of the particle before it is scattered. For 
higher MF a particle will be more quickly 
accelerated and it can attain a higher energy 
before it may escape to the far-upstream 
region of the forward shock. Therefore, 
the MF determines the maximum energy 
of both protons and electrons, the latter in 
addition suffering energy losses determined 
mainly by the same MF. 

The actual mechanism of MF amplifica- 
tion (MFA) is a subject of ongoing research. 
Despite recent progress 6^, 2a, 35 1, it is not 
yet clear what mechanism (or most prob- 
ably a combination of mechanisms) is re- 
sponsible for induction of MHD turbulence 
on scales needed for the particle accelera- 
tion (for a recent review see 



61 and refer- 



ences therein). Due to the uncertainty in 
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Figure 2: Magnetic field profiles in the SNR for D 
and P type models at different ages. 



the dominant mechanism of MFA, and the 
difficulty of taking into account the com- 
plexiety of the MFA process, the magnetic 
field in SNRs is often phenomenologically 
parametrized. If the MP is assumed to be 
frozen in, the energy density of the field may 
correlate with the thermal and relativistic 
particle energy densities. Here we investi- 
gate the impact of widely used MP proffies 
on the spectra of accelerated particles. We 
introduce a set of simple scalings which de- 
scribe the temporal and spatial evolution of 
the magnetic field in a SNR. The models are 
listed in Tabled] and are explained below in 
subsection 14.21 

Prompted by recent high resolution X- 
ray observations, which claim to detect non- 
thermal emission from the reverse-shock re- 



Table 1: Nomenclature of MF models used in this 
study. 



Profile 

MF ^^"^-^ 


Density 


Pressure 


75/iG 

-|- Alfvenic drift 


D075I 
D075A 


P075I 
P075A 


300/iG 

-|- Alfvenic drift 


D300I 
D300A 


P300I 
P300A 



gion of some SNRs [38|, |39 



we have 

chosen to consider particle acceleration at 
the reverse shock in addition to that at the 
forward shock. The major argument against 
the ability of the RS to accelerate parti- 
cles is that an exceptionally low MF exists 
in the ejecta, a theoretical argument ob- 
tained under the assumption that the mag- 
netic flux at the surface of progenitor star is 
conserved during the expansion of the SNR. 
However, the detection of non-thermal X- 
rays from the RS implies that there may be 
some processes acting against rapid dilution 
of the MF, and that a large MF at the re- 
verse shock can exist. These processes may 
be the fluid-dynamical MHD turbulence in 
the ejecta and CD region - the initial ac- 
celeration of CRs, while the MF was high 
enough, with subsequent streaming through 
a confined ejecta region (accelerated CRs do 
not leave the upstream region of the RS as 
easy as that of the FS). It is also not clear 
whether the MF in the ejecta, to be scaled 
according to flux conservation, must be the 
fleld at the surface of SN progenitor. Fi- 
nally, the fact that shocks in SNRs are col- 
lisionless implies the existence of a flnite MF 
which mediates the collisionless shock. Ad- 
ditional justiflcation for the presence of MF 
at the RS can be found in 44 . 
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4-2. Parametrizations of Magnetic Field 
Profiles 

The actual amplification of the mag- 
netic field is a complicated, turbulent pro- 
cess that is not well understood. Many 
published studies of particle acceleration 
in SNRs assume that the MF is propor- 
tional to the density of the plasma (e.g., 
33 . 44 . 62 . 52II). An alternative approach 



is generally adopted in models of SNe radio 
emission, where the MF energy density in 
the remnant is proportional to the thermal 
energy density, and thereby the gas pres- 
sure |63| (for review see [64]). Following 
these scalings, we set up 8 MF models, 4 'D' 
models following the density profiles and 4 
'P' models scaling with the square root of 
pressure profiles. All our models are time 
dependent, but the matching at the FS is 
different for the two classes of MF models. 
D models assume a constant value of the 
MF at the FS, Bq = Bps, and only the 
scaling inside the SNR proportionally fol- 
lows the time-dependent density distribu- 
tion, B{r,t) = BQp{r,t)/ p{Rps,t). In con- 
trast, P models assume that the MF at the 
FS and inside the SNR both evolve with 
time according to the square root of pres- 
sure distribution, B{r,t) = a/ C87rP(r, t). 
Here, p{r,t) and P{r,t) are the density and 
the pressure profiles of the plasma inside 
the SNR, and C is a proportionality con- 
stant chosen so that B^Rpsytm) = Bps, 
where tm = 500 years is the mid-time of 
simulations. We assume two possible val- 
ues of the MF at the SNR forward shock. 
Bps = 75 fiG and Bps = 300 pG. These 
values are taken in accordance with aver- 
age limits given by the observational data 
on the MF for young SNRs in X-rays 471 
and gamma-rays 46[. The proportionality 



for Bps = 300 /iG. The amplified value 
of the MF at the RS varies between ^ 5 
and ~ 30 pG depending on the age and 
model. Even if we start from a rather di- 
luted MF, this requires amplification factors 
of ~ 50 — 300, which are not unreasonable 
in CR-invoked or fluid-mechanical turbulent 
amplification. 

The hydrodynamic profiles allow for MF 
scaling only inside the shocked region. Gen- 
erally one can expect that the MF amplifi- 
cation processes operate throughout the en- 
tire shock precursor region j26|, and thus 
the MF should show a smooth transition 
from far-upstream values to amplified val- 
ues, which are then compressed to the given 
amplitudes at the shock [Bps in D models 
and B{Rps,t) in P models). We therefore 
assume that the MF falls off to the inter- 
stellar field (5 fiG) at the distance of 5% of 
the FS radius ahead of the FS, and down 
to a very small ejecta field (0.01-0.1 /iG) 
at the distance of 5% of the RS radius to- 

This transition may 



ward the interior 44 



constant, C, is rather small in both cases, 
C = 0.19% for Bps = 75 pG and C = 3% 



be explained by an exponential decrease in 
the number density of high energy particles, 
which may cause MF amplification. In gen- 
eral, one can also expect that the diffusion 
coefficient becomes larger than Bohm diffu- 
sion, albeit still much lower than the Galac- 
tic one due to the enhanced flux of CRs. 
The transition to the Galactic diffusion oc- 
curs at the escape boundary, where the con- 
centration of CRs tends to zero. Here we as- 
sume that Bohm diffusion operates up to es- 
cape boundary, Resc- In our simulations, we 
assume that in young SNRs the MF is fully 
isotropic on account of efficient turbulent 
field amplification, therefore the MF com- 
pression ratio is aMF = -\/0^H-"2a2)73 = 
a/TT, where a is the gas compression ratio. 
The distributions of the MF for different 



model types and different times are given 
at Fig El 

4-3. Alfven Drift of Scattering Centers 

As noted earlier, we also investigate the 
influence on the particle spectra of a possi- 
ble drift of the scattering centers upstream 
of the shocks. Resonant streaming insta- 
bilities in the precursor region of the shock 
generate outward-moving Alfven waves 261 
that may have a significant phase velocity 
in the presence of amplified MF. As par- 
ticle are scattered on excited MHD waves, 
the Alfven velocity of the waves may effec- 
tively change the compression ratio felt by 
the particles, and the spectrum becomes no- 



ticeably softer [65|, l32|, |35|, |33| . In cases in 



volving an Alfvenic Mach number Ma ^ 10, 
the test-particle slope is softer, s < 2, than 



the classical solution [3J]. It should be 
noted that non-resonant instabilities often 
involve turbulence with very small phase ve- 
locity [25|] for which Alfvenic drift is ignor- 
able. Which type of instability dominates 
and what the turbulence properties are in 
the nonlinear regime is the subject of ongo- 



ing research [66|, |67|, |68| . Here we paramet 



rically study the impact of Alfvenic drift. 
In half of our models we introduce Alfvenic 
drift in the upstream regions of both shocks 
assuming that the effective velocity of the 
scattering centers equals the Alfven veloc- 
ity, va = Bps/ V^'^PiSM, plus the upstream 
gas velocity. In the downstream region, 
we assume that the scattering centers are 
isotropized and there is no Alfvenic drift. 
We mark the models including Alfvenic drift 
with "A " while the others are marked with 

5. Results and Discussion 

In this section we present our main re- 
sults. We discuss spectra of particles ac- 



celerated at the reverse and the forward 
shocks, analyze their properties and evolu- 
tionary differences, and explore the intrin- 
sic differences in particle spectra imposed 
by different parametrizations of the MF and 
how the Alfvenic drift of scattering centers 
affects particle spectra. Finally, we calcu- 
late the emission from the particles in SNR 
and discuss the resulting observable proper- 
ties which could be compared to the exper- 
imental data. 

5.1. Particle Spectra 

5.1.1. Forward vs. Reverse Sfiock 

The non-thermal spectra of accelerated 
particles are characterized by a power-law 
index and a cut-off energy. The power- 
law slope according to DSA theory de- 
pends solely on the compression ratio of 
the plasma flow (velocity jump across the 
shock). At the location of the shock, where 
the plasma is compressed, particles adiabat- 
ically gain energy and then lose it every- 
where in the expanding flow. In the evolv- 
ing system, the average compression ratio 
should be considered. Our hydrodynamical 
simulations were performed with sufficiently 
high resolution that the deviations of the 
shock compression ratio do not affect the 
spectral slope through the entire evolution. 
The compression ratios in our simulation 
are aps = 4.01±0.09 and aRs = 3.97±0.11, 
for the FS and the RS respectively. 

In the models without Alfvenic drift, 
the particles accelerated at the FS obey 
the classical simple power-law distributions 
with index s ~ 2 up to the energy, at which 
the escape of the particles from the system 
(for protons) or the energy losses (for elec- 
trons) inhibit further acceleration. As pre- 
dicted by DSA [69|, for protons Emax,p ~ 
V^f^Bot. Therefore, the increase of Emax,p 
due to system expansion (higher energies 
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Figure 3: The FS (thick lines) vs. the RS (thin hnes) contributions to the spectra of protons (left) and 
electrons (right) for different SNR ages. 



are needed to leave the system) is eventu- 
ally compensated by the decrease of Vsh, so 
Emax,p decreases after reaching a maximum. 
This is illustrated by the time evolution of 
FS proton spectra in Fig. [3] (left), which 
shows the D300I model, in which Bq is kept 
constant, thus isolating Emax,p oc V^f^t. For 
electrons, Emax,e evolves in the same man- 
ner up to the time when the energy gain 
per cycle no longer exceeds the energy loss 
per cycle due to synchrotron radiation, i.e, 
tacc = trad- A secoud, typically lower, char- 
acteristic energy, Erad, is approximately set 
by equality of the radiative-loss timescale 
with the age of the remnant. Above E^adi 
the electron spectra assume a quasi steady- 
state, changing slowly on account of the 
evolution of the gas-flow profiles. Besides, 
above Erad the electron spectra steepen by 
one power approximately (Fig. |3l right). 

The ability of the RS to accelerate par- 
ticles depends on the MF configuration at 
the location of the RS. With the MF profiles 
assumed here and discussed above, it is ob- 
served that the RS is capable of accelerating 
particles to sufficiently high energies and in- 
tensities to contribute to the total volume- 



integrated particle spectra. These results 
are at least consistent with many of the ob- 
servations noted above. In order to obtain 
the volume-integrated spectra, N{E), one 
performs an integral of the spatially differ- 
ential particle distribution, N{x,E), over 
the entire volume of the simulation grid up 
to Resc = '^Rfs'- 

N{E)=A'kR%s I N{x,E)x^dx (6) 
Jo 

where x = r/Rps- In all plots featur- 
ing particle spectra we show the quantity 
N = N{E)/At[R],s. In the absence of 
Alfvenic drift, the RS spectra show slopes 
similar to those at the FS, with the ex- 
ception of some hardening at higher ener- 
gies. The maximum energies for protons are 
lower than at the FS on account of a lower 
MF strength, which allows particles to es- 
cape the acceleration region at lower ener- 
gies. In contrast, the maximum energies for 
electrons are nearly as high as at the FS 
mainly because a weaker MF also reduces 
the rate of radiative losses of electrons. The 
intensity of the volume-integrated RS spec- 
tra falls with time relative to the FS spectra 



as the number of injected particles becomes 
smaller. Also the ratio of the RS surface 
to FS surface decreases with time, which 
means fewer particles that could be poten- 
tially injected cross the RS. 

The stronger spectral hardening and 
more significant bump at high energies 
make the RS volume-integrated spectra dif- 
ferent from those of the FS. The spectral 
hardening is most visible at early evolution- 
ary epochs and gradually disappears with 
time for both shocks, but at a somewhat 
faster rate for the FS. The bump in the for- 
ward and the reverse shock spectra is lo- 
cated close to Emax, which is beyond Erad 
for electrons. The presence of these spec- 
tral features and their evolution can be ex- 
plained by the geometry of our spherically- 
symmetric system. 

In order to get volume-integrated spec- 
tra we integrate over the entire simulation 
domain, so both upstream and downstream 
regions contribute to the volume-integrated 
spectra. The volume-integrated upstream 
spectra usually have hard slopes, s ~ 1, 
and the intensity of the upstream spectra 
around Emax is comparable with that of 
the downstream spectra. Thus, in spectra 
integrated over both upstream and down- 
stream regions one can see the contribution 
from the hard upstream spectra in high- 
energies close to Emax, whcrcas the con- 
tribution of low-energy particles is negligi- 
ble. Besides, the i^max-region of the up- 
stream spectra exhibits a bump which ap- 
pears on account of the exponential de- 
crease of the MF in the upstream direc- 
tion. As a result, the diffusion time of par- 
ticles located close to the shock is longer 
than that of particles further away from 
the shock. Therefore, the recently accel- 
erated particles around E^ax are accumu- 
lated near the shock in the upstream re- 



gion. The spectral hardening and the bump 
is strongly visible in the RS spectra (Fig. [3l 
right) because the volume of the upstream 
regions of both shocks are significantly dif- 
ferent in the spherically-symmetric geome- 
try assumed here, contrary to the case of 
plane-parallel shocks. This also explains the 
different time evolution of the spectral fea- 
tures. The upstream region of the FS ex- 
pands significantly faster than the upstream 
region of the RS, thus the CR number den- 
sity in the upstream region of the FS de- 
creases faster and contributes less to the 
total volume-integrated spectra. We note 
the presence of a peculiar wiggle in the elec- 
tron spectra. The MF in the upstream re- 
gion is lower than downstream, resulting in 
weaker synchrotron losses for the electrons 
and therefore Erad — Emax,e- Consequently, 
the upstream electron spectra have a sig- 
nificantly higher intensity near Emax,e than 
the downstream electron spectra. There- 
fore, in the total volume-integrated spectra 
the bump stands out beyond Erad in the 
Emax,e rcgiou. For protons, Emax,p in the 
upstream and downstream spectra nearly 
coincide, and a small bump appears near 
Emax,p in the total volume-integrated pro- 
ton spectra. Again, due to the faster de- 
crease in the CR number density ahead of 
the FS with time, the feature is more promi- 
nent in the RS spectra. 

5.1.2. Influence of Magnetic Field Profiles 
and Strengths 
Most important for DSA is the configura- 
tion, strength, and evolution of the MF at 
the location of the shocks and their precur- 
sor regions. Since here we adopted similar 
profiles of MF in the precursor regions of 
the shocks, only a different time evolution 
of the MF strength at the location of the 
shocks may change the shape of the volume- 
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Figure 4: The time evolution of the spatial distribution of protons (top row) and electrons (bottom row) at 
energy of 20 TeV. 



integrated spectra. It was noted, however, 
that the differences between models in the 
spectra of particles created by the RS and 
the FS are marginal, in fact only Erad and 
Emax truly depend on the MF strengths. 
While for protons a higher MF would mean 
a higher Emax,p, for electrons this would 
mean a lower Erad and Emax,e due to ra- 
diative losses. 

The MF profiles in the SNR interior affect 
the radial distribution of CR particles and 
their subsequent emission, especially so for 
leptons. In FigJH we plot the evolution of 
the radial distribution of protons and elec- 
trons at the energy of 20 TeV, which could 
be a characteristic energy of the particles 
producing high-energy emission, discussed 
in more detail in subsection 15.2.21 We use 
models without Alfvenic drift as it does 
not change the distribution significantly - 
it would only slightly reduce the intensity 



at the FS. One can clearly see that elec- 
trons suffer significant losses in the shocked 
region where the MF is high. The num- 
ber density of protons falls nearly exponen- 
tially from the FS towards the CD. There is 
a noticeable peak in the CR distributions 
near the RS. Whereas electrons peak ex- 
actly at the RS, protons are accumulated 
between the RS and the CD. At 100 years, 
the CR number density close to the RS is 
even higher than at the FS. At 400 years 
the peak is still clearly visible, however al- 
ready at 1000 years only a minor bump is 
observed. In the ejecta region, where the 
MF is low, both particle species show nearly 
uniform distributions. Although the D075I 
model follows the trends described above, it 
is a clear outlier. The cut-off energy of par- 
ticles accelerated by the RS in this model 
is around 4 TeV, so one finds at 20 TeV 
almost no particles from the RS. This ex- 
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Figure 5: The spectra of the FS (thick hnes) and 
the RS (thin hnes) for different SNR ages in P300A 
(top) and D300A (bottom) models. 

plains the absence of a peak near the RS and 
the low CR number density in the interior 
of the SNR. Besides, the electrons have suf- 
fered little energy losses in the initial phases 
of SNR evolution {t < trad) compared with 
other models (remember that P075I model 
have initially a high MF). 

5.1.3. Impact of Alfvenic Drift 

We have explored the influence of 
Alfvenic drift on the particle spectra of the 
SNR. Interesting is that it affects our mod- 
els in a different way during the SNR evo- 
lution. We observe that the FS and the RS 
spectra in P models are noticeably softer 
during all evolutionary times. This is be- 
cause in P models the Alfvenic Mach num- 
ber, M^, is almost constant for both shocks 
during the evolution. The MF strength, 
-B, evolves as the square root of the pres- 



sure, -^/P, and hence roughly scales with 
the shock velocity. Therefore, spectra show 
nearly the same softening at all times as 
shown at the top of Fig. |5l Contrary to P 
models, in D models the effect of Alfvenic 
drift is different for the FS and the RS. At 
the FS B is kept constant during the evolu- 
tion, and Ma is very large at early times be- 
cause the Alfven velocity, va-, is small com- 
pared to the shock speed. At later times 
Va remains constant, but the shock speed 
drops, and so Ma becomes small. This 
implies that the softening of the FS spec- 
tra in D models was small at early times 
and increased during the evolution of the 
SNR. In contrast to the FS, B is evolving 
at the RS as density, p. The initially large 
eject a density caused a high Alfven velocity 
{va ~ ^/p in D models) and rather small 
M^, though larger than at the FS, which 
caused some softening of the spectra. Over 
time p dropped significantly, and therefore 
Ma became large. Therefore, the softening 
of the RS spectra in D models is stronger at 
earlier times and almost not noticeable at 
later times (Fig. El bottom). Obviously, all 
these effects and characteristics are stronger 
in models with high MF strength. 

5.1.4. Total spectra from both shocks 

The shape of the total volume-integrated 
spectra depends on the SNR age on account 
of the changing importance of the contribu- 
tions from the forward and reverse shocks. 
For very young Type la SNRs, and keep- 
ing in mind the MF profiles assumed herein, 
a significant contribution to the total spec- 
trum arises from the RS. However, already 
after a few hundred years the bulk of par- 
ticles is provided by the FS. Nevertheless, 
particles accelerated at the RS continue to 
create pecuhar spectral features up to late 
times of the SNR free expansion phase (see 
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Fig. |6]). Added to this is the different im- 
pact of Alfvenic drift on the spectra pro- 
duced with different MF models and the 
time evolution of this impact, especially in 
D models. The total spectra are very dif- 
ferent from a single-shock test-particle solu- 
tions as seen at Fig. [61 We plotted 6 out of 
8 models because D075A and P075A mod- 
els show similar shapes as D075I and P075I 
models, apart from being slighter softer. 

At 100 years the contribution of the RS to 
the total spectra is significant in all models. 
With age this contribution decreases. Af- 
ter 400 and 1000 years, the proton spectra 
in P models look like broken power laws or 
vary gradually, much smoother than expo- 
nentially cut-off. Since the RS proton spec- 
tra in D models have a somewhat smaller 
cut-off energy and are harder than the FS 
spectra, if Alfvenic drift is assumed, the RS 
maintains a significant contribution to the 
total spectrum until late epochs. Because of 
it and because of the bump contributed by 
the upstream region of the FS, the proton 
spectra of D models show some small-scale 
concavity at high energy through all times, 
however at later times it is shallower and 
broader (bumps in both RS and FS spectra 
become less significant). 

The total electron spectra (Fig. |H1 bot- 
tom row) also have interesting features. On 
account of energy losses, the contribution of 
the RS in the cut-off region is stronger. This 
creates either small-scale concavity or a bro- 
ken power-low in the loss-affected region 
(between Erad and Emax,e)- The stronger 
the contribution from the RS is, the more 
pronounced a concavity is observed. When 
the contribution of the RS is mild, then one 
sees a broken power-law. In some models, 
(i.e., D300I and D300A), the shallow con- 
cavity in the loss-affected region is visible 
throughout all late times. It is striking that 



the electron spectra in single-shock NDSA 



simulations [70|] have similar concavity in 
the loss-affected spectral band. Addition- 
ally, one can note that near E^ad the spectra 
of P models are smoother than the spectra 
of D models. 

We must make an important remark con- 
cerning particle spectra. The plotted spec- 
tra are integrated over the entire simulation 
domain. However, one finds (15.2. ip that the 
main contribution to the emission is done 
by the particles confined inside the remnant. 
These particle spectra are typically softer at 
high energies because they do not include 
the hard upstream particles. For illustra- 
tion see the thin lines of D075I models in 
Fig. El which show separately SNR-confined 
(downstream) particles and peaked at high 
energies upstream particles. Only inverse 
Compton emission because of low MF in 
the upstream region can benefit from the 
upstream particle population. 

5.2. Particle Radiation 

We consider three basic non-thermal 
emission processes, synchrotron and inverse 
Compton (IC) emission of electrons, and 
neutral pion decays originating in collisions 
of CR protons with protons at rest. The 
synchrotron emission is calculated accord- 
ing to 7l|] . We consider only the CMB pho- 



ton field for the IC scattering. Therefore 
we may use a simple rescaling of the syn - 
chrotron emission to IC scattering 72|, [73 



We calculate the gamma-ray emission pro- 
duced by pion decays according to the pro- 
cedure described in 7^. Thermal emission 
is not considered. However, we made es- 
timates that even assuming instant equili- 
bration the maximum of the thermal emis- 
sion does not exceed 10% of the synchrotron 
emission at the respective energy for all con- 
sidered ages. In particular, at 3 keV the 
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Figure 6: The time evolution of volume-integrated proton (top row) and electron (bottom row) spectra for 
different models and SNR ages. Thin lines for the D075I model in the top row display the downstream and 
upstream components of the proton spectra. 



thermal component may comprise ~ 3% of 
the X-ray emission of the models considered 
here. The reader should be aware though 
that we used the same injection efficiency 
for electrons and protons which results in 
high amplitudes of leptonic emission. If the 
electron injection efficiency is lowered signif- 
icantly, below 3 keV the brightness profiles 
may be affected. 

5.2.1. Radiation Spectra 

The radiation spectra show diversity dur- 
ing SNR evolution (see Fig. [7]). In general, 
the radiation spectra are much smoother 
than the parent particle distributions. The 
leptonic spectra show broad features, which 
are the result of energy losses experienced 
by electrons during their acceleration his- 
tory. A noticeable feature of the leptonic 
spectra is concavity, which appears in dif- 



ferent models at different times. At 100 
years it is strongly visible in P300I and 
P300A models, but in D075I and D075A it 
is nearly invisible. Later, at 400 years, it 
appears in D300I and D300A models while 
it disappears in D075I and D075A, and it 
starts to vanish in P300I and P300A mod- 
els. At 1000 years it is visible only in D300I 
and D300A models. The concavity reflects 
the change in electron slope on account of 
the contribution of the RS between E^ad 
and Emax,e. It interesting to see that the 
volume- integrated synchrotron and IC spec- 
tra have different shapes. While the syn- 
chrotron emission comes mainly from the 
shocked region, a signiflcant fraction of IC 
emission comes from the ejecta and up- 
stream region with low MF, where the elec- 
tron spectrum is different. The pion-decay 
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Figure 7: The time evolution of volume- integrated emission from the SNR due to synchrotron (top row), 
inverse Compton (middle row) and pion-decay (bottom row) for different models and SNR ages. 



spectra are much smoother than their lep- 
tonic counterparts and can be characterized 
by a power-law with a very gradual and 
broad cut off at relatively low energy be- 
tween a few 100 GeV to a few TeV. This 
makes the spectra rather soft in the en- 
ergy band observed with modern Cherenkov 
telescopes. The pion-decay emission comes 
mainly from inside the SNR, with a small 
amount of upstream contribution. 

5.2.2. Brightness Profiles 

To permit morphological studies of SNR 
at different ages, we computed radial inten- 



sity profiles at characteristic wavebands for 
current instrumentation (radio @1.4 GHz, 
X-rays @3 keV, gamma-rays @1 TeV). In 
Fig. |H]we separately plot the emission of lep- 
tonic and hadronic origin. The profiles are 
normalized to the emission maximum along 
the SNR radius. The absolute intensity may 
be obtained from the spectra given at Fig. [71 
We plot here "I" models only. Their "A" 
counterparts show similar profiles, but the 
intensity of the FS is a bit lower, making 
the profiles look somewhat flatter. 

It is astonishing to see how the brightness 
profiles undergo severe changes with time. 




Figure 8: The radial intensity profiles of the SNR in different models at different wavelength and ages. Radio 
@1.4 GHz (dotted-pink), X-rays @3 keV (solid red), IC @1 TeV (long-dashed green), pion-decay @1 TeV 
(dashed blue). 



At 100 years both shocks are bright and 
visible in all bands except for the FS in ra- 
dio waves, which are emitted by low-energy 
electrons. The contribution of the RS to the 
total emission is high at the initial stages. 
The low-energy electrons cannot propagate 
far from their point of origin, and therefore 
the relatively high injection rate at the RS 
in the early phase renders it much brighter 
than the FS in radio band. Pion-decay 
gamma rays are copiously produced in the 
ejecta region on account of the high target 
density there. At 400 years the contribution 
of the RS region becomes less prominent, 
but is still visible. Only in synchrotron X- 



rays (and IC in D075I/A models) the RS be- 
comes insignificant. The 3 keV synchrotron 
X-rays are created by electrons of rather 
high energy (~ 13.5 TeV for B = 300 fiG 
and ~ 27 TeV for B = 75 /xG). At this 
age the RS is no longer able to boost elec- 
trons up to these energies. The electrons 
accelerated at the FS suffer from losses on 
their way to the SNR interior, and there- 
fore the synchrotron emission becomes no- 
ticeably filamentary, especially in high-MF 
models. A characteristic feature develops 
in IC emission of the models with high MF. 
A plateau of high intensity appears in the 
ejecta region and a low intensity between 
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the two shocks. In the shocked region, high- 
energy electrons are few on account of se- 
vere energy losses in the high MF, so the 
IC emission is suppressed. At the same 
time, the electrons that propagated into the 
ejecta region, where the MF is low, accumu- 
late there and account for bright IC emis- 
sion. This trend continues with age, and af- 
ter 1000 years we observe a bright plateau of 
IC emission from the SNR center, while the 
FS appears very filamentary. There is al- 
most no trace of emission from the RS after 
1000 years in the other wavebands. X-ray 
synchrotron emission is seen as a thin fila- 
ment near the FS, especially in the models 
with high MF. 



6. Conclusions 

We have studied particle acceleration by 
both forward and reverse shocks in young 
type-la Supernova Remnants. Rather than 
starting with Sedov models, as is most com- 
monly done, we used gas- flow profiles in 
the ejecta-dominated phase, derived from 
1-D hydrodynamical simulations, to solve 
the cosmic-ray transport equation in a test- 
particle approach. We analyzed how differ- 
ent phenomenological parametrizations for 
the magnetic-field distribution and strength 
affect the maximum energy and spatial dis- 
tributions of accelerated particles. Addi- 
tionally, we have explored the influence of 
possible Alfvenic drift of scattering centers 
on the total particle spectra in the upstream 
region of both SNR shocks and find our re- 
sults for the forward shock in agreement 

We in- 



with previous findings [32|, |35|, 134 , 
vestigated the properties of the resulting 
non-thermal radiation, i.e., its spectra and 
the spatial distribution of the intensity at 
characteristic wavelengths. 



We have demonstrated that it is impor- 
tant to account for the contribution of the 
RS to cosmic-ray particle population in the 
initial 400 - 600 years of SNR evolution, 
given the widely used parametrizations of 
the MF that are assumed here. At that 
time the RS is able to accelerate particles 
up to very high energies, and the number of 
accelerated particles is comparable to that 
at the FS. This is well visible in volume- 
integrated particle spectra, the distribution 
of high-energy particles in the SNR, as well 
as in the morphology of particle radiation. 
The significance of the RS contribution falls 
with time. Although the RS contribution is 
still visible in particle spectra up to 1000 
years or so, it is barely noticeable in the 
emission spectra and morphology, which is 
in agreement with observations. 

We found that the total volume- 
integrated particle spectra of the SNR 
can be very different from single-shock 
test-particle solutions, especially for very 
young SNRs. The choice of MF profiles 
affected the relative contribution to the 
total spectrum of the RS and the FS, and 
consequently the emission spectra in the 
high-energy region varied. Additionally, 
Alfvenic drift of scattering centers in the 
shock precursors gives rise to additional 
spectral feature, whose appearance and 
time evolution depends on the choice of 
MF profile. Generally, P models are more 
sensitive to Alfvenic drift than D models. 
Obviously, the models with high MF are 
affected stronger than those with low MF, 
for which the effect is subtle. 

The particle distributions computed here 
show a variety of spectral shapes. Interest- 
ing is that some of the spectra (both elec- 
tron and protons) exhibit features which 
may be considered similiar to ones shown 
using NDSA, such as spectral concavity and 
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high-energy bumps. However, concavity in 
our spectra arises in a much smaller high- 
energy band compared to the broadband 
concavity of NDSA. Besides, it is not re- 
flected in pion-decay spectra as these are 
produced predominantly by confined par- 
ticles in the SNR. Additionally, our test- 
particle spectra are softer at high energies 
than those of NDSA, even if both are mod- 
ified by Alfvenic drift, which may provide 
better agreement with gamma-ray spectra 
observed from SNRs. 
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